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G^ ! Abstract 

We have investigated numerically the quantum evolution of a (5-like 
C^ I wave-packet in a quenched disordered medium described by a tight- 

binding Hamiltonian with long-range hopping (band random matrix 
i-A . approach). We have obtained clean data for the scaling properties in 

P3 I time and in the bandwidth b of the packet width M and its fluctua- 

^ ' tions Ajy,j with respect to disorder realizations. We confirm that the 

fluctuations of the packet width in the steady-state show an anoma- 
lous scaling A^^/M ~ b~^ with 6 = 0.75 it 0.03. This can be related 
kS I to the presence of non-Gaussian tails in the distribution of M. Fi- 

C^ • nally, we have analysed the steady state probability profile and we 

have found 1/6 corrections with respect to the theoretical formula de- 
rived by Zhirov in the b ^ oo limit, except at the origin, where the 
corrections are 0(l/v6). 

PACS numbers: 05.45.Mt, 71.23.An, 72.15.Rn, 05.60. Gg 



1 Introduction 

Band random matrices (BRM) represent an effective model for both ID dis- 
ordered systems with long-range hopping and quasi-lD wires[||]. The band- 
width b plays the role of the range of the interaction in the first case, the one 
of the square root of the number of independent conduction channels in the 
second. Up to now, studies have been mostly devoted to the analysis of the 
stationary solutions of the Schrodinger equation and to the corresponding 
spectral properties of BRM's |0]. Much less is known about the solutions of 
the time dependent Schrodinger equation, a topic on which only a few studies 
have been performed |^, ^. The partial analogy of this latter problem with 
the 'dynamical localization' phenomenon in the kicked rotor [^] suggests that 
an initial delta-like packet spreads diffusively and eventually saturates to a 
localized state. The width of this asymptotic packet for BRM's is of the 
order of 6^ lattice sites, i.e. the same order as the localization length of all 
the eigenf unctions [0. 

The theoretically predicted scaling laws for the mean square displacement 
M were tested numerically and a comparison of the asymptotic form of the 
wave-packet with a theoretical formula [^ derived for the ID Anderson model 
was attempted [Q. More recently some new theoretical results appeared 
which give a formula for the time asymptotic packet in the BRM model in 
the large b limit 0. Therefore, it became important to check numerically this 
formula and to both investigate how the packet reaches its time asymptotic 
shape and measure the size of finite b corrections. For what the time evolution 
is concerned some phenomenological expressions were suggested in ref. [Q, 
based on a power-law convergence of the mean square displacement to its 
steady state value. However, the presence of large statistical fluctuations 
prevented the authors of ref. [Q from assessing whether the time asymptotic 
scaling is ruled by power law corrections or by the logarithmic corrections 
to the t~^ dependence suggested by rigorous results obtained for the ID 
Anderson model 0. The fluctuations A^q- of the width of the asymptotic 
wave-packet with the realization of the disorder constitute an even more 
controversial issue, since not even the scaling behaviour is clearly understood. 
Some evidence of an anomalous behaviour was presented in two previous 
studies of the same problem ||, ^ and in the kicked rotor p. In all cases 
the numerics was too poor to make a convincing statement about the value 
of the anomalous exponent. 

The bottleneck of the previous simulations was the slowness of the inte- 



gration scheme, a 4-th order Runge-Kutta with a small time step to obtain 
a good conservation of probability over a long time span. This low efficiency 
prevented from reaching sufficiently large values of b and from considering 
a large enough number of realizations of the BRM's. We have instead im- 
plemented a 2-nd order Cayley algorithm, which, being unitary, exactly con- 
serves probability, although the one-step integration error is larger than the 
one of the Runge-Kutta scheme (a situation similar to those of symplectic 
algorithms in classical Hamiltonian dynamics). This has allowed us to more 
than double the maximum bandwidth (from b = 12 to 6 = 30) and to increase 
the statistics by a factor four (in the worst case). 

As a result, we have been able to complete an accurate analysis of the time 
evolution of the mean square displacement, finding that there is no need to 
invoke effective formulas with a power-law time dependence even at relatively 
short times. We have found a clean evidence of an anomalous scaling of the 
relative fluctuations of the packet width, which behave as Aj^/M ~ b~^ with 
6 = 0.75 ± 0.03. In order to con&m this anomaly, we have investigated the 
statistics of the packet width at a specific time in the localization regime. The 
probability distributions at various b values, when appropriately rescaled, 
superpose, and the resulting universal (6 independent) curve is definitely 
different from a Gaussian with an exponential tail at large M values. Finally, 
we have compared our results with the theoretical formula for the asymptotic 
wave-packet finding a convincing agreement. The finite b corrections to 
the b ^ oo Zhirov expression are of order (l/b). 

2 Model and numerical technique 

We have considered the time-dependent Schrodinger equation 

'-gr= 2^ Hiji^j (1) 

j=i-b 

where ipi is the probability amplitude at site i and the tight-binding Hamil- 
tonian Hij is a real symmetric band random matrix. The band structure of 
the Hamiltonian is determined by the condition 

Hij = if |z - jl > 6, 



the parameter b setting the band-width; the matrix elements inside the band 
are independent Gaussian random variables with 

{Hi,)d = and {{H,,f)d = l + 5^j 

where the symbol {■)d stands for the average over different realizations of 
the disorder. In the present work we have considered the evolution of an 
'electron' initially localized at the centre (identified with the site i = 0) of 
an infinite lattice. To this aim we have analysed the solution of equation (|lD 
corresponding to the initial condition 

^. (t = 0) = 6io. 

Since the wave-packet evolves in a supposedly infinite lattice, it is neces- 
sary to avoid any spurious boundary effect due to the inevitably finite size of 
the vectors used in the numerical computations. This goal has been achieved 
by resorting to a self-expanding lattice, i.e. a lattice whose size is progres- 
sively enlarged according to the development of the wave-function. At each 
integration step, our program checks the probability that the electron is in 
the leftmost and rightmost b sites, adding 106 new sites whenever the ampli- 
tude \ipi\ is larger than e = 10"'^ in at least one of the 26 outermost sites. We 
have separately verified that e is small enough not to significantly affect the 
computation of the probability distribution. For instance, by lowering e by 
an order of magnitude, the mean squared displacement (computed over the 
same disorder realizations) changes only by a few percent. Since this system- 
atic error is not larger than the uncertainty due to statistical fluctuations, it 
is not convenient to reduce the cut-off as it would turn out in a slower code 
with a consequent reduction of the statistics. 

The Schrodinger equation (|I]) was integrated by approximating the evo- 
lution operator exp {—iHt) with the Cayley form 

, .,^,, l-iH5t/2 

which implies that the values of the wave-function at two successive time- 
steps are related by 

1 + \iH5t\ i){t + 5t) = (l - \^H5t\ i){t). (3) 



Solving the band diagonal system of equations (|^) allows one to determine 
ipit + St) once ip{t) is known. Cayley's algorithm is a standard tool for the 
computation of the solutions of the Schrodinger equation (see for instance 
ref. [0); to the best of our knowledge, this is the first application to the 
specific field of random Hamiltonians with long-range hopping. Cayley's 
form (^ for the evolution operator has two relevant features: it is second- 
order accurate in time and unitary; in addition, the corresponding integration 
scheme (0) is stable. Stability is essential in order to study the long time 
evolution of the wave-packet; as for unitarity, it ensures the conservation of 
probability and, together with second-order accuracy in time, allows one to 
choose time steps 6t two or three order of magnitude bigger than those used 
in Runge-Kutta integration schemes. Indeed, we could make use of a time 
step 6t ~ 10~^, to be compared with the time step 6t ~ 10~^ — 10~^ used for 
the same problem in Refs. 0, ^ . To ascertain how large a 6t could be used, we 
have compared the solutions obtained through Cayley's algorithm at various 
6t with the exact solution of the Schrodinger equation (P, computed by 
diagonalizing the Hamiltonian (to avoid boundary effects due to the finite size 
of the diagonalized matrices, we have considered sufficiently short evolution 
times). By this way we came to the somewhat surprising conclusion that the 
validity range of the approximate equality (@) extended up to time steps as 
big as 6t ~ l/v6 (the scaling of 5t with the band-width b was necessary to 
compensate the opposite scaling of the energy eigenvalues with v6)- To check 
this conclusion, we have computed the mean squared displacement in the 
localized regime for several values of 5t in the range 10~^/v6 — 1/Vb, finding 
differences of a few percent, not larger than the statistical fluctuations. 

This can depend on the fact that the long time evolution of the wave- 
packet seems to be led by the eigenstates at the band centre. Indeed, in the 
energy representation, the exact evolution operator and the Cayley form can 
be written as 



exp {-iH6t) = Y, \n)e-'^-^' 

n 



n\ 



l-ziJ5t/2_^,^,_,^„(,,) 



1 + iH5t/2 

with 

(/)„(5t) = 2 arctan (E„(5t/2) , 

where \n) is the eigenvector corresponding to the energy E^- These equa- 
tions show that, for increasing 5t, the approximate equality (^ holds true 



only in the subspace spanned by the eigenvectors corresponding to the band 
centre, while at the band edge, where the eigenvalues \En\ tend to y/b, the 
coefficients exp {—iEnSt) and exp {—icpn^t) become quickly different. There- 
fore, the eigenstates at the band edges appear to play a minor role in the 
time evolution. This is probably due to the shorter localization length of 
such states compared with those in the centre: in fact, for the mean square 



displacement, all eigenstates are weighted with their localization length |11 



3 Results 

To investigate the time evolution of the wave-packet, we have computed the 
mean square displacement 

M{b,t) = {u{t))^/ £ fmt)A ■ (4) 

\j=-oo / ^ 

Previous studies of this problem strongly suggest that M satisfies the 
scaling relation 

M{b,t)=b^M{T = t/b^/^) , (5) 

for large enough values of the bandwidth b. Nevertheless, in Ref . [Q , where 
the most detailed numerical investigation has been carried out, it was not 
possible to obtain a clear verification of the scaling law (||) due to the poor 
statistics and the small values of b. 

The faster integration algorithm described in the previous section has 
allowed us both to average over more realizations (400 in the worst case), 
i.e to reduce statistical fluctuations, and to reach larger values of b (namely, 
b = 30 instead of 6 = 12 as in Ref. IQ). The results reported in Fig. |l] for 
several values of b are clean enough to show a convincing convergence from 
above to a limit shape. In other words, there is no possibility to interpret 
the deviations as a signature of a different scaling behaviour for M. 

In order to perform a more quantitative analysis, we have proceeded in the 
following way: M{t, 6)Qhas been averaged over the time interval 20 < r < 30 
to obtain the more statistically reliable quantity {M)t{b). By assuming a 
dependence of the type 

(M),(6)=M^(l + ar"), (6) 



^We have added the variable b to underhne the residual but asymptotically irrelevant 
dependence on the band-with. 
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Figure 1: Rescaled mean squared displacement M vs. rescaled time r 
t/h^''^. From top to bottom: 6 = 8, 12, 16, 22, 26, 30 



we have fitted the three parameters Moo, a and a, finding that the conver- 
gence rate a is very close to 1 (0.95), i.e. that the finite-band corrections 
are of the order 1/6. The fitted value of M^o is 0.61. The results for the 
finite-band correction 



SM = M^- {M)t{h) 



(7) 



are plotted in Fig. 0. 

The good quality of our numerical data suggests also the possibility to 
compare the temporal behaviour with the available theoretical formulas. In 
particular, it has been argued in Ref. that the existence of the so-called 
Mott states should imply a (lnt)/t convergence of M to its asymptotic value. 
Therefore, we propose the following expression 



M{T,h) = M{oo,b) 1- 



l + Alnjl + T/tpY 
1 + r/tD 



(8) 
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Figure 2: Convergence towards zero of the finite band correction (|^). The 
slope is 0.95 thus close to -1. 



which is the simplest formula that we have found able to reproduce also the 
initially linear (i.e. diffusive) regime. For each value of b, the best fit is so 
close to the numerical data of Fig.^ to be almost indistinguishable from them 
(this is why we do not report the fits on the same figure). The meaningfulness 
of the above expression is further strengthened by the stability of the three 
free parameters M{oo,b), A and tc, which allows the calculation of a b- 
independent diffusion constant (see below). From the values oi M{oo,b), we 
can extrapolate the asymptotic value M(oo, cxd) in exactly the same way as 
we have done for (M)t(6), finding once more a (l/6)-convergence to a value 
around 0.70. This result is to be compared with the theoretical prediction 
M{oo, oo) ^ 0.668 DTSf. The deviation of about 0.03 can be attributed to the 
accuracy of the integration algorithm. 

Another important parameter that can be extracted from formula 



IS 



the diffusion coefficient. Indeed, by expanding Eq. (^ for small r, we find 

MM = MlESii!) (1 _ ^) = c (9) 

r to 

The diffusion constant D turns out to be close to 0.50 for all values of b and, 
what is more important, close to the value that we obtain from a quadratic 
fit of the initial growth rate of the packet. This is a very encouraging result, 
since it confirms the correctness of formula (§) for both the diffusive and the 
localized regimes. Let us notice that the value D ~ 0.50 is somewhat smaller 
than the one reported in |^ (D ~ 0.83). Taking into account statistical 
fluctuations and systematic deviations, we find that D = 0.50 ± 0.05. 

In past papers, a phenomenological expression involving a power-law con- 
vergence to the asymptotic value of the mean squared displacement has been 



proposed P, 12|, arguing that it should provide an effective description of 



both the diffusive and localized regime 

M{t, h) = M(oo, b) (l - \ ] . (10) 

The success of expression (H) shows that there is no need to introduce anoma- 
lous power laws to reproduce the numerical findings. However, for the sake 
of completeness, we have fitted our numerical data also with Eq. (|T0|), finding 
an equally good agreement. Therefore, on the basis of the quality of the fit 
we cannot conclude which of the two expressions is better; nevertheless, it 
is worth recalling that the former one has the correct asymptotic behaviour 
and, moreover, the fitted parameters are more stable. 

A much more controversial situation exists about the fluctuations of the 
packet width. Let us introduce the r.m.s. deviation 

Am(&, t) = ^Hty), - {u{t))l . (11) 

In fact, it has not yet been clarified how the above variable scales in the 
large-6 limit. In particular, the correct value of the scaling exponent u in the 
relation 

A,^(6,t) = 6^AM(r = t/6='/^6) (12) 

is still unknown; this is why the b dependence in Am is explicitly maintained. 
A scaling like b^ would imply that the packet-width is not a self-averaging 
quantity, since the relative size of the fluctuations would not go to zero 

9 



for increasing b. Conversely, an exponent u = 3 corresponds both to self- 
averaging and 'normal' behaviour. In fact, the number N^ of independent 
channels (lattice sites) actively contributing to the localized region (i.e. the 
localization length) is of the order 6^. If we assume that all such contributions 
to the second moment M are independent of one another, then we are led to 
conclude that the relative fluctuations should decrease as l/y/TT^ = 1/6, thus 
yielding an absolute growth as b^. Since previous studies @] have suggested 
a small anomaly, i.e. u slightly larger than 3, we have chosen to report 
the behaviour of Am(t, 6) for u = 3. The data shown in Fig. ^ reveals a 
drastically different behaviour from what observed in Fig. [1|. First of all, the 
curves tend to grow for increasing b; moreover, there is no obvious indication 
of a convergence to some finite value. Altogether, these features imply that u 
is strictly larger than 3, qualitatively in agreement with previous simulations. 
In order to perform a more quantitative analysis, we have computed the 
average of Aj^ and rescaled it to the average of M, 

A^(6) = ^^ . (13) 

^^ ' {M)t 

{{■)t is again to be interpreted as the average over the time interval 20 < r < 
30). The advantage of this renormalization, already adopted in Ref. Q], is 
that it reduces finite-band corrections. The results reported in Fig. |, reveal 
a clean power law decay with an exponent 6 ~ 0.75. This value is slightly 
larger than the one found in the previous studies, but follows from a much 
cleaner numerics. A more global check of the scaling behaviour can be made 
by plotting the rescaled fluctuations 

Al{T)^b'^. (14) 

for the various values of b. The optimal value of 6 can thus be identified as 
that one yielding the best data collapse. The curves reported in the inset of 
Fig. ^ have been obtained for 6 = 0.75. It is necessary to modify 6 by at 
least ±0.03 units in order to see a significant worsening of the data collapse. 
Accordingly, the best estimate of the anomalous exponent is 5 = 0.75 ±0.03, 
so that the dependence of Am on b in Eq. (^) is removed for u = A— 6 ~ 3.25. 
In order to find further support for this anomalous behaviour, we have 
investigated the probability distribution P{M) for the second moment at 
the time r = 30 (the longest time we reached for the larger 6- values), i.e. 

10 
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Figure 3: Rescaled fluctuations of tlie mean squared displacement Am in 
formula (|12[) with z/ = 3 vs. rescaled time r = tj}?!'^ . The values of 6 
increase from bottom to top 6 = 8, 12, 22, 26, 30. There is a clear tendency 
to grow for increasing h values. 



when the wave-function has almost entered the steady-state regime. The 
construction of reliable histograms has forced us to consider smaller values 
of h. In fact, we have studied the cases 6 = 8, 6 = 12 and h = 22, using 
10^ realizations of the disorder in the first two cases and 10^ in the last one. 
The results are reported in Fig. |^, where, following a method suggested in 
Ref. |T^, we have conveniently rescaled the probability distribution P[M). 
In particular, defining by Mav the average value of M at r = 30 and a{M) the 
standard deviation over the ensemble of disorder realizations, we have plotted 
P'{M') = a{M)P{M) vs. M' = {M - Mav)l(y. After this rescaling, the 
distributions P'{M'), corresponding to the three b values, have zero average 
and unit standard deviation. It is remarkable to notice that all curves nicely 
overlap indicating a striking scaling behaviour. A further important feature 
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Figure 4: Averaged and normalized fluctuations of the mean squared dis- 
placement in formula (p!3| ) vs. b (log- log scale). The line is a fit with a power 

(see Eq. (^ for its 
0.75. 



law b with 6 = 0.75. In the inset we report A^(r 



definition) for the same values of b as in Fig. ^ and 6 



is the deviation from a Gaussian behaviour, especially for large values of M', 
where a clear exponential tail is visible. The dotted line just above the three 
curves (corresponding to the pure exponential exp(— M')) has been added to 
give an idea of the decay rate which is slightly larger than 1. The results of 
this analysis are important in two respects: i) the exponential tail 'explains' 
the difficulties encountered in getting rid of statistical fluctuations in the 
estimate of Am] ii) the deviations from a Gaussian behaviour provide an 
independent evidence of the anomalous scaling behaviour of the fluctuations. 
It is interesting to remark that a preliminary quantitative comparison has 
revealed a striking identity of the probability P'{M') with the distributions 
found in Ref. [|14| for such diverse quantities as the magnetization in the 2D 
XY model and the power consumption in a confined turbulent fiow. This 
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close correspondence deserves further investigations. 
.0 




Figure 5: Scaled probability distribution of M' = (M — Mav)/cr. The dotted 
line has been drawn to guide the eye to the exponential (non- Gaussian) 
behavior. 



Finally, we want to compare our results with the theoretical predictions 
for the asymptotic shape of the wave-packet. In |^ a reasonable agreement 
was found between the numerical data and the formula obtained by Gogolin 
for strictly one-dimensional systems p. Since a theoretical expression has 
been derived in the meantime for quasi-lD systems 0], it is desirable to 
compare our data also with this expression. In Fig. ^ we present the disorder- 
averaged probability profiles {\4'j(t)\'^)d = f{j,t) for large times, rescaled 
under the assumption H] 



fs{x) = b'^f{j,oo); 



X 



j/b' 



(15) 



and compare them with Zhirov's theoretical formula, which is denoted by 
the white line. No appreciable deviation is noticeable except for the extreme 
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part of the tails, where it is reasonable to expect numerical errors due to 
boundary effects. 
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Figure 6: Probability profile rescaled using formula (^) for several b values 
vs. X = j/y^, compared with Zhirov's theoretical prediction (white line). 



The good overlap is partly due to the (unavoidable) choice of logarithmic 
scales in Fig. |^. However, if we zoom the region around the maximum (with 
the exception of the zero channel) , one can see in Fig. |^ a slow tendency of the 
various curves to grow towards the theoretical expectation. This is consistent 
with the behaviour of M reported in Fig. |I|, which reveals a convergence 
from above for the mean squared displacement. It is interesting to notice 
that all such deviations are mainly due to the finiteness of b while the lack 
of asymptoticity in t appears to be much less relevant. 

Finally, we consider separately the zero channel, i.e. the return probabil- 
ity to the origin /s(0). In Fig. ^ we plot this quantity versus r for different 
values of b. In all cases, a quite fast convergence, as compared with the 
behaviour of the packet-width, to the asymptotic value is clearly seen. In 
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Figure 7: Zoom of the central region of Fig. ^. Zhirov's formula is now the 
black line. 



practice, as soon as r is about 1, the average value of /s(0) reaches the 
asymptotic value. It is instructive to compare our numerical findings with 
the asymptotic (in time and b) analytic expression /^(O) = 6. By fitting 
the dependence of the time average of fg (in the interval 1 < r < 30) on 
b as in Eq. ^, we find that the asymptotic value is about 5.7 and that the 
convergence rate is 1/Vb. The numerical value is in a reasonable agreement 
with the theoretical one, considering that it is the result of an extrapolation 
of data already affected by errors of the order of a few percent. The non 
trivial part of the result is the rate of convergence of this probability, which 
is definitely slower than the 1/6 behaviour displayed by the second (and other 
low order) moments. The behaviour of the return probability, however, is 
in agreement with the theoretical predictions made in 0, where the finite 
b deviations from the asymptotic steady-state probability distribution were 
estimated to be of order 0(1/a/6) in the |a;| < 1/6 neighbourhood of the 
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origin. 




100 



Figure 8: Return probability to the origin /s(0) vs. r for several b values: 
from bottom to top 6 = 8, 16, 22, 30. 



The 6-dependence of the return probability is further illustrated in Fig. P, 
where we plot the deviation 5/ = 6 — /s(0) from the asymptotic value 
limb^oo/s(0) = 6 as a function of b in bilogarithmic scale. The displayed 
numerical values were obtained by averaging the return probability both 
over disorder realizations and the time interval 20 < r < 30; the data were 
then fitted with two expressions, exhibiting deviations from the asymptotic 
value /s(0) = 6 of order 0(1/6) and 0(1/6°^) respectively. In the second 
case, the exponent a was used as a fitting parameter and the best fit value 
was a = 0.53. As can be seen from Fig. ^, the power law with 0(l/v6) 
corrections fits the data much better than the one with deviations of order 
0(1/6). 
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Figure 9: Deviations of the return probability from the asymtpotic value vs. 
b (log- log scale). The circles represent numerical data, the dashed line is a 
fit with a power law 1/6, the continuous line is a fit with a power law 1/6", 
with a = 0.53. 



4 Conclusions and perspectives 

We have studied the time evolution of an initial 5-like wave-packet in a ID 
disordered lattice with long-range hopping. The main results of this pa- 
per are the following. We have confirmed with clean numerics the scaling 
law (H) for the mean square displacement M, first proposed and studied in 
Ref. [^]. This scaling law is valid in the large b limit; here we have found 
that finite b corrections are of the order 1/6. We have proposed formula 
(H) for fitting the time evolution of M towards its steady state value; this 
formula contains the logarithmic corrections suggested by the existence of 
Mott states. We confirm the presence of an anomaly in the scaling law of the 
relative fluctuations Aj^^/M of the mean square displacement, finding that 
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they vanish for large b as 6~°-^^. We have hnked this anomaly to the presence 
of non- Gaussian fluctuations of the mean square displacement. In fact, the 
probability distribution of M displays an exponential tail for large values of 
M. The conveniently rescaled probability strikingly coincides with the distri- 
butions obtained in Ref. ||14| for such diverse quantities as the magnetization 
in the 2D XY model and the power consumption in a confined turbulent 
flow. The degree of universality of such distribution deserves further investi- 
gations. Finally, we have compared the numerical results on the steady state 
probability profile with the theoretical formula proposed by Zhirov for large 
b, finding a good agreement. We have computed for the first time finite b 
corrections, obtaining 0(1/6) deviations for the moments of the probability 
profile and 0{1/Vb) corrections for the return probability to the origin. 
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